Numerical treatment of long-range Coulomb potential with Berggren bases 
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The Schrodinger equation incorporating the long-range Coulomb potential takes the form of a 
Fredholm equation whose kernel is singular on its diagonal when represented by a basis bearing 
a continuum of states, such as in a Fourier-Bessel transform. Several methods have been devised 
to tackle this difficulty, from simply removing the infinite-range of the Coulomb potential with 
a screening or cut function to using discretizing schemes which take advantage of the integrable 
character of Coulomb kernel singularities. However, they have never been tested in the context 
of Berggren bases, which allow many-body nuclear wave functions to be expanded, with halo or 
resonant properties within a shell model framework. It is thus the object of this paper to test 
different discretization schemes of the Coulomb potential kernel in the framework of complex-energy 
nuclear physics. For that, the Berggren basis expansion of proton states pertaining to the sd-shell 
arising in the A ~ 20 region, being typically resonant, will be effected. Apart from standard 
frameworks involving a cut function or analytical integration of singularities, a new method will be 
presented, which replaces diagonal singularities by finite off-diagonal terms. It will be shown that 
this methodology surpasses in precision the two former techniques. 

PACS numbers: 02.60.Nm, 02.70.Hm, 23.50.+Z 
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I. INTRODUCTION 

The infinite-range of the Coulomb potential has always 
been a source of numerous problems in both theoreti- 
cal and numerical studies of the quantum mechanics of 
charged particles. 

From a theoretical point of view, the 1/r Coulomb 
asymptotics implies the presence of an essential singu- 
larity of the Coulomb Green's function at k — in the 
complex plane, contrary to that of potentials quickly van- 
ishing atr-> +00 [l| . This has prevented for a long time 
the use of the simple and elegant demonstration of the 
completeness of bound and scattering states solutions of 
the one-body Schrodinger equation of R. Newton as 
the latter explicitly demands r 2 \V(r)\ to be integrable 
on [0 : +00 [, where V{r) is the basis-generating potential. 
Until recently, no simple proof existed for the Coulomb 
case, as the only available methods requested abstract 
Lebesgue measure theory Q , thus hindering any physical 
understanding of this fundamental property. Hence, new 
methodologies had to be devised in order to have both 
simple and rigorous demonstrations, based on a general- 
ization of the method of R. Newton [J] or on closeness ar- 
guments with standard complete sets of states Q , which 
rely heavily on the analytical properties of Coulomb wave 
functions. 

At the numerical level, the computational effort neces- 
sary to evaluate Coulomb wave functions is tremendous 
compared to that of its uncharged counterpart, the ana- 
lytical Ricatti-Bessel function. Indeed, even though ana- 
lytically expressible in terms of confluent hypergeometric 
functions |6j, Coulomb wave functions are very difficult 
to calculate, especially in the complex plane @, Hj|- They 
can vary by many orders of magnitude over a relatively 
small region and thus demand the use of many different 
numerical techniques @, Q. Added to that, they exhibit 



a cut in the complex plane on ] — 00 : 0] , which creates a 
supplementary difficulty as the functions issued from con- 
tinued fractions, for example, are analytical therein @,@]- 
When one considers a Fourier-Bessel transform of poten- 
tials bearing a Coulomb tail, no problem appears at the 
theoretical level. This comes from the fact that the latter 
can always be decomposed as Vo(r) + c/r, where Vo(r) is 
a quickly vanishing potential and c is a constant. Indeed, 
on the one hand, the Fourier-Bessel transform of Vo(r) is 
well-behaved, and on the other hand, the term propor- 
tional to 1 jv has a Fourier-Bessel transform analytical for 
all orbital momentum i, equal to Q^[(k 2 + k' )/(2kk')]/n, 
where Qi(x) is the Legendre function of the second kind 
&. However, Qi(x) diverges like ln(l — x) when x — > 1 
9], so that the momentum space representation of V(r) 
has a logarithmic singularity on its diagonal at k = k', 
and consequently cannot be discretized readily, as the 
diagonal of the discretized matrix is infinite. Neverthe- 
less, the Fourier-Bessel transform of 1/r is integrable, so 
that it is possible to solve this problem with subtraction 
frameworks, where singularities are integrated analyti- 
cally, leaving a regular rest integral to be handled nu- 
merically [HI • This has been successfully applied for the 
diagonalization of the Coulomb potential in momentum 
space in Refe. 0,113, E2- 

However, no attempt has been made using Berggren 
bases [l3[ . The latter are sets of wave functions generated 
by finite-range potentials and comprise bound, resonant, 
and scattering states of complex-energy, which are solu- 
tions of the one-body Schrodinger equation. Berggren 
sets of eigenstates are widely used in nuclear physics 
to solve the many-body problem, as they allow com- 
plex nuclear wave functions to be expanded in a ba- 
sis of Slater determinants built from Berggren one-body 
states, in the so-called Gamow shell model framework 
(see Ref.[l4| for a recent review on that subject). Due 
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to the presence of scattering states in Berggren bases, 
the asymptotic behavior of complex many-body wave 
functions can be precisely reconstructed, which is im- 
possible to attain in practice with bases of well-bound 
states, such as sets of harmonic oscillator states. They 
are then very well suited to expand halo and unbound 
many-body nuclear states [15rl25f , whose properties can- 
not be understood without a precise reproduction of their 
asymptotic properties. Generating potentials are usually 
of Woods-Saxon type, which mimic the effect of the in- 
ert core on valence nucleons (2(| . Hartree-Fock potentials 
are also utilized, as their variational character allows con- 
figuration mixing to be minimized [l9l [22| ■ In all these 
cases, wave functions cannot be conveyed in an analyt- 
ical form and have to be computed numerically. Com- 
pleteness properties of Berggren sets of states have been 
established firstly by T. Berggren for uncharged parti- 
cles [HI, while the charged case could be handled only 
recently in Ref.@. Completeness relations have been 
tested by numerically diagonalizing Hamiltonians sus- 
taining spherical and axially-deformed nuclear potentials 
for both charged and uncharged particles [H, [ill H3, [H| , 
but where the Coulomb potential is absent from the diag- 
onalized kernel. If one considers nuclear wave functions 
of a fixed nucleus, this problem is just anecdotic, as, due 
to its one-body character, the infinite-range part of the 
Coulomb potential, which is spherical, can always be in- 
cluded in the basis-generating potential, leaving a finite- 
range residual interaction to diagonalize [25j . However, 
this is no longer possible if one wants to calculate ob- 
servables connecting nuclei bearing different numbers of 
protons, as the same Berggren basis has to be used to ex- 
pand the many-body wave functions of both father and 
daughter nuclei. Problems also occur to calculate isospin 
operators expectation values with Berggren bases (see 
for example Ref . (25| for a study of isospin mixture of iso- 
baric analog states of A = 6 nuclei). Indeed, it is needed 
therein to evaluate one-body matrix elements (a\t \b), 
where \a) and \b) can be scattering wave functions. If pro- 
ton and neutron wave functions of a given partial wave 
are generated by the same potential, (a^jfe) matrix ele- 
ments are proportional to a Dirac delta distribution. In 
the opposite situation, however, Dirac delta distributions 
cannot arise as proton and neutron radial wave functions 
of same quantum numbers (£, j) are not orthogonal, so 
that matrix elements are undefined. Consequently, both 
proton and neutron Berggren set of states have to be gen- 
erated by the same one-body Hamiltonian, for example 
with the basis-generating potential of neutron states, so 
that the infinite-range Coulomb Hamiltonian has to be 
diagonalized along with the residual nuclear interaction. 

Thus, in this paper, a one-body Hamiltonian possess- 
ing an infinite-range proton potential, thus presenting a 
Coulomb point potential asymptotic, will be diagonal- 
ized with Berggren bases employing three different dis- 
cretization schemes. The first scheme is simply to cut 
the Coulomb potential at a finite distance R, determined 
so as to provide optimal results. The second scheme is a 



subtraction technique, akin to that expressed in Ref. [121] . 
and the third scheme consists of replacing the diverging 
diagonal elements by off-diagonal terms, which become 
diagonal at the continuum limit. These methods will 
be described in detail in Sec. [Hi where they will be re- 
spectively labeled by "the cut method," the "subtraction 
method," and the "off-diagonal method." Numerical ap- 
plications will be presented in Sec. Mil for three different 
partial waves, S1/2, ^3/2 and d 5 / 2 , where both narrow 
and broad resonant states will be considered. The con- 
sidered expanded one-body wave functions, namely the 
lsi/2, 0g?3/2 and 0e? 5 / 2 proton states, arise typically for 
nuclei for which A <~ 20, so that their study is impor- 
tant for future calculations of many-body states of light 
nuclei with the Gamow shell model. Numerical methods 
proper to the use of proton Berggren bases states, gener- 
ated by a potential with a Coulomb asymptotic, will also 
be discussed. The conclusion is stated in Sec. IIVI 



II. THEORETICAL BACKGROUND 

A. Expression of the Hamiltonian kernel in a 
Berggren basis 

One considers a Berggren basis of a partial wave of 
quantum numbers whose completeness relation 

reads Hal: 



u n (r)u n (r')+ / u k {r)u k {r') dk = 6(r-r'), (1) 



nG(h,d) 



where L + is a complex contour of momenta k starting 
at k = and ending at fc = +00, representing the scat- 
tering part of the Berggren basis constituted by the \uk) 
states, and where the discrete sum runs over bound (b) 
and resonant, or decaying (d) states \u n ), the latter res- 
onant states being situated between the real fc-axis and 
L + (see Ref.[l8[ for details). The energy of the discrete 
states \u n ) will be denoted as e n and that of scattering 
states as e^. Discrete states \u n ) are normalized to 
one and scattering states \uk) to Dirac delta [18j. They 
are generated by a one-body Hamiltonian of the form: 
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h=^-+V ws {r) + V c {Z c ,r) 
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Vws(r) = -V f(r) - 4(1*. s)Ko : 
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m = 

V c (Z c ,r 



1 + cxp 



r-Ro 



C c Z c erf(ar) 



(2) 

(3) 
(4) 



where m is the effective mass of the proton, Vws{t) is 
a potential of Woods-Saxon type, with its diffuseness d, 
radius Ro, central and spin-orbit depths V and V so , re- 
spectively, I is the orbital momentum operator and s the 
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spin operator, and V c is the Coulomb potential, propor- 
tional to an error function, defined with the Coulomb 
constant C c , the charge acting on the proton Z c and a 
constant a. The radial function entering the Coulomb 
potential V c (Z c ,r) is standard and arises from the use 
of a Gaussian charge density for the closed core [29M3ll j . 
For simplicity, the finite-range part of both basis and di- 
agonalized Hamiltonians will be taken to be the same so 
that they will differ only through their charge Z c . We 
will denote the charge of the basis potential as Z^ , that 
of the diagonalized potential Zc, and their difference, 
which will enter the kernel to diagonalize, will be denom- 
inated as AZ C = Zc d ^ — Zc b \ The eigenstate \<j>) of the 
diagonalized Hamiltonian is to be expanded with Eq. ([1} : 

|0) = E c n \u n ) + / c k \u k ) dk, (5) 

ne(b,d) L+ 

where the coefficients c n and c k have to be determined. 
Its energy will be designated as E. Using Eqs. (JI]) and 
([2]) , the Fredholm equations which must be solved to as- 
certain the c„ and c k coefficients are readily inferred: 

c n e„ + ^ c n >(u n >\V c (AZ c ,r)\u n ) 

n'e(b.d) 

+ / c k > (u k >\V c (AZ c ,r)\u n ) dk' 
Jl+ 

= E Cn Vn € (b, d), 
c k e k + E c n >{u n ,\V c {AZ c ,r)\u k ) 

n'e(b.d) 

+ / c k <{u k ,\V c {AZ c ,r)\u k ) dk 1 
Jl+ 

= E c fe Vfc e L+. (6) 

Matrix elements of Eq. © have to be calculated using a 
complex rotation of r [U [HJ due to the unbound char- 
acter of basis states: 

(u a \V c {AZ c ,r)\u b ) = [ u a (r)V c (AZ c ,r)u b (r) dr 
Jo 

+ E el9 / <«(z(x))-^—^u" b »(z(x)) dx, (7) 
e —± J R + xe w 

where \u a ) and \ub) are two Berggren basis states, R is 
a radius after which complex rotation is effected, suf- 
ficiently large to sustain erf(or) ~ 1 and /(r) ~ 
in Eqs. (j3]) and (j4|), ui = ± characterizes the asymp- 
totic character of wave functions components, as u(r) = 
u+{r) + u-{r) = C+H+(kr) + C-H- ri {kr) for r > R, C± 

being a constant and H^(z) a Coulomb wave function 
of orbital momentum £ and Sommcrfcld parameter 77, of 
outgoing (+) or incoming (— ) character, 9 is the rota- 
tion angle, which depends on w a and u b and is chosen 
so that the associated improper integral converges and 



z(x) — R + xe l6 . Coulomb wave functions are imple- 
mented using recently developed techniques mixing ana- 
lytical formulas and direct integration t 8j. Complex ro- 
tation of r is also employed to normalize the discrete 
states \u n ) [see Eq. ([IJ] [18|. One has to pay attention if 
k(R + xe 10 ) crosses the negative real axis with complex 
rotation of r, as H^(z) becomes discontinuous therein. 
Indeed, Coulomb wave functions bear a cut by defini- 
tion, whereas the wave functions u^(r) are everywhere 
continuous, as solutions of a differential equation of sec- 
ond order. Crossing cannot occur in practice for reso- 
nant states, sustaining sufficiently small width-to-energy 
ratios, while non-crossing can be enforced for scatter- 
ing states by taking a L + contour sufficiently close to 
the real A:-axis. It is, however, unavoidable for bound 
states when cos(6>) < 0. In the latter case, a linear 
combination of H^ n (kr) and H^ ! (kr), whose coefficients 
are simple functions of rj and k [8|, has to be added to 
the initial wave function u(r) = u + (r) = C + H^ ri (kr) 
to suppress discontinuity. Because, for bound states, 
\H7 (k(R + xe l9 ))\ -!• when x +00 with cos(6>) < 0, 
this modification does not change the asymptotic prop- 
erties of u(r). 

Due to analyticity of u a (r) and u b (r) in the complex 
r-space, the matrix element embedded in Eq. ([7]) is in- 
dependent of R and 9, so that it is a genuine defini- 
tion of operators acting on Berggren basis states. How- 
ever, as one might expect, Eq. ([7J cannot be used for 
\u a ) = \ub) = \u k ). Although improper integrals for 
which uj a uj b = 1 arc always well defined with complex 
rotation of r, no rotation angle 9 can be found to have 
improper integrals converged when Lu a uj b = —1. This is 
straightforward to demonstrate using the asymptotical 
expressions of u^(z) for |z| — > +00, which reads: 

u ± (z) ~ exp[±i(kz - r}ln(2kz))], (8) 

up to an unimportant constant. As a consequence, sim- 
ilarly to the representation of 1/r in momentum space, 
diagonal matrix elements are infinite (i.e. they cannot be 
regularized with complex rotation of r) and off-diagonal 
matrix elements (u k \V c (AZ c , r)\u k >) for which k 7^ k' , 
where complex rotation is always available, diverge like 
ln(fc — k') when k' — > k. Modification of diagonal matrix 
elements in Eq. ([7]) is then necessary to discretizc Eq. ^ . 



B. Discretization of Berggren basis 

The most efficient method to discretize Eq. © is to 
utilize a Gauss-Legendre quadrature for the L + contour 
of Eq. © [ljj. In order to obtain a symmetric matrix in 
the resulting eigenproblem, discretized scattering states 
are set equal to yfwl u ki (r), where is the weight as- 
sociated to the ki abscissa. The discretized Berggren 
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completeness relation then reads: 
^2 u n {r)u n (r') 

n£(6,d) 
Not-l 

+ ^ u k z {r)u kl (r') Wi ~ 5(r - r') 

i=0 

iV r e S -l 

<=> Ui(r)ui(r') 

i=0 
N-l 



+ 51 Ufci = ( r ) u ki s ( r ') w ls ~ S(r - r') 

i=N rcs 
N-l 

^ ^ Ui(r)ui(r') ~ S(r - r'), 



(9) 



where N res is the number of bound or resonant states, 
Nql is the number of Gauss-Legendre points, N — 
N res + Nql is the total number of basis states, i s = 
i — N res , notation which we will maintain from now on, 
and Ui(r) is a bound or resonant state u n (r) if i = n € 
[0 : N res — 1] and the discretized state y/Wi a Uk ia (r) if 
i £ [N res : N — 1] . The energy of the i-th state will be 
designated as e^. The resulting discretized eigenstate im- 
plemented from diagonalization of discretized kernel will 
be naturally denoted as: 



N-l 
i=0 



Cut method 



(10) 



This approximation is the crudest, as it simply removes 
all improper integrals in Eq. (JT)), which is equivalent to 
replacing V c (Z c ,r) of Eq. Q by V c {Z c ,r) He(R - r), 
where He is the Heaviside function. As a consequence, 
Eq. ([7]) becomes: 



(u a \V c (AZ c ,r)\u b ) 



u a (r)V c (AZ c ,r)u b (r) dr (11) 



This removes all singularities in the Hamiltonian kernel of 
Eq. ([2]) , but at the price of introducing a cut-dependence 
on R which has to be assessed. Its discretized eigenprob- 
lem is then straightforward to convey utilizing Eq. (|10[) 
(see Sec. Ill Bl for notations): 



N-l 

Ci ei + 2J Ci>(ui<\V c (AZ c , r)\ui) = E c l; 

i'=0 

< i < N. 



D. Subtraction method 



(12) 



This is the standard method to treat integrable sin- 
gularities in Frcdholm kernels II Oil . The main idea is to 



separate the kernel integral into two parts, one whose 
integrand is regular and the other singular but analyt- 
icall y in tegrable. However, contrary to the situation of 
Ref. [12lj . it is necessary to terminate the L + contour at 
finite k = k max . Indeed, allowing the L + contour to go to 
+oo in the fc-plane, with Gauss-Laguerre quadrature for 
example, would demand the consideration of wave func- 
tions of very high linear momenta, consequently tremen- 
dously oscillating. The latter would have to be used in 
radial integrals which have to be implemented numeri- 
cally, such as the nonanalytical kernels of Eq. © , hence 
inducing numerical instability. Because, in practice, con- 
vergence with k max is rather fast, it is prefcrrable to uti- 
lize Gauss-Legendre quadrature on a L + contour ending 
at k — k ma x ■ 

In order to apply the subtraction method, the integral 
of Eq. ^ involving k and k' is rewritten the following 
way: 



Ck'(u k '\Vc(AZ c ,r)\u k ) dk! 



L + 



(c k > Uk'(r) V c (AZ Cl r) u k (r) 



Ck s k '{r) Cc s fe (r) I dr 



dk 1 





C c AZ C 


/> 






r 



Sk ) dk', 



(13) 



where \sk) is a sine function [i.e. (r\sk) — \/2/7rsin(fcr)], 
so that one has just added and subtracted the I = 
Fourier-Bessel transform of the Coulomb potential 
(C c AZ c )/r. The calculation of the radial integrals in- 
side the first fc'-integral of the right-hand side of Eq. (fTS"]) 
can always be performed with complex rotation of r 
when k ^ k' (see Sec. Ill A[) . When k' ~ k, the con- 
sideration of the asymptotics of u^(r) for r — > +oo [see 
Eq. flSJ)], the fact that 2nC + C~ — 1 (arising from Dirac 
delta normalization of scattering states [18[) and that 
Ck' - c k ~ (fc' - k)(dck"/dk")k"=k (cfe is an analytic 
function of k) imply that the radial integrals problem- 
atic with complex rotation of r; that is those for which 
k = k' and Lu a u>b = —1 (see Sec. Ill Al for notations), are 
convergent in Eq. (|13[) . 

The second integral of the right-hand side of Eq. (Ti~3"]) 
can be evaluated analytically: 



Sk' 

L+ 

C c AZ C 



a az c 



Sk ) dk' 



[\n(k + k') - ln(fe - k')] dk' 



k 

C c AZ C 
n 



[ln(fc + k') - \n{k' - k)\ dk' 



[(k 

in ax 

+ k) \n(k max + k) 



[kmax k)ln{k max k) 2k\n(k)] 



(14) 
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where the fact that the complex linear momentum fc be- 
longs to the L + contour has been taken into account. The 
discretization of Eq. ([6]) can now be effected, which re- 
sults in a different treatment of diagonal and off-diagonal 
terms involving scattering states (see Sec. Ill Bl for nota- 
tions): 



zoidal rule, one obtains: 

1X fn|fc-fc'| dk' 
+ I " ln(fc' - k) dk' 



ln(fc - k') dk' 



Ik 

N-l 



M 



Ci 6z ~t~ 



N-l 

E 

i'=0 



+00 



Ci'iui'lVdAZ^r^Ui) = E a , < i < N re 



< J [ Ui(r) 2 V c {AZ Cl r) 
C AZ 

Ci \(kmax ~\~ ) hl(fcj, 



sfe ls (r) s 



C, AZ ( 



dr 



ln(fc - iAk) Afc + ^ ln(iAfc) Afc 

i=0 i=l 

Afc Afc 

— ln(fc) — \n(k rnax - k) 

[N ln(Afc) + ln(iV!) + M ln(Afc) + ln(M!)] Afc 

Afc Afc 

m ( fc ) 7T \n(k max - k) 



2 
(fc, 



2 

fc) ln(fc r , 



7T 



{J^max fci, 

C c AZ C 



1 ln(fc„ 



- Ca 



a AZ C 



i s -l 

E 

i'=0 



Ngl-1 

E 



ix - k is ) - 2k is ln(fc; 8 )] 
[ln(fc is +fci/)-ln(fc ia -h',)] 

Wi> [ln(fc ia + h> ) - ln(fcj> - fc ls )] 



+ fcln(fc) - k - Afc In 
Afc -> 0. 



- fc) - (fc 
Afc 
2^T 



0(Afc 2 ), 



(16) 



If one replaces the infinite value of In | fc — fc' | at fc' = fc 
by ln[Afc/(2?r)] in the first line of Eq. (fill), its end 
point contribution of the trapezoidal rule of the inte- 
grals defined on [0 : fc] and [fc : k max ] at fc' = fc is 

k„ 



N-l 

E 

i'=0 



equal to (Afc/2)ln[Afc/(27r)], so that / ln|fc-fc'|c2fc' 
(15)_ 7o 

is reproduced up to Afc with the trapezoidal rule, as 
Afcln[Afc/(27r)] is then added to Eq. (fig]). As a con- 
sequence, we introduce the transformation fc — > fc — 
Afc/(47r) and fc' -> fc + Afc/(4?r) at fc = fc'. This has 
the advantage not to spoil the precision of the inte- 
gral if one adds a symmetric function /(fc, fc') regular at 
fc = fc' to In |fc - fc'|, as /(fc - Afc/(47r), fc + Afc/(47r)) = 
/(fc, fc) + 0(Afc 2 ), immediate from the second-order Tay- 
lor expansion of /(fc, fc') at fc = fc'. Note that 0(Afc ) is 
the error one expects from the standard trapezoidal rule. 

Even though one cannot assess analytically the error 
made using Gauss-Legendre quadrature instead of the 
trapezoidal rule, it will be proved to be very small nu- 
merically. For that, we consider the following integral: 



E. The off-diagonal method 



/(fc, fcmaa?) 



(s k ,\V c {AZ C) r)\s k ) dk', (17) 



The exchange of the diagonal infinite matrix elements 
by close off-diagonal matrix elements is motivated by the 
analytical approximation of the integral of In \k — k'\ be- 
tween fc' = and fc' = k max acquired from the trape- 
zoidal rule. For that, one uses two grids of equally space 
points on [0 : fc] and [fc : k max ], of same discretization 
step Afc (k m ax is assumed to be an integral multiple of 
Afc), so that N = fc/Afc and M = (k max — fc)/Afc are 
their respective number of points. Removing the infinite 
value at fc' = fc from the sums deduced from the trape- 



which is the integral entering Eq. ([5]) possessing the same 
singularities as that of Eq. (|16[) . where basis functions 
Uk (r) have been replaced by sine functions s k (r) and the 
complex L + contour is hereby the real segment [0 : k max \. 
If one takes V c (AZ c ,r) to be exactly equal to the point 
Coulomb potential (C c AZ c )/r [i.e. a — > +00 in Eq. @], 
Eq. (JTTJ) is analytical [see Eq. flU}]. When a is finite, 
it is nevertheless possible to evaluate almost exactly the 
integral of Eq. (fTT]) rewriting V C (AZ C , r) as [V C (AZ C , r) - 
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(C c AZ c )/r] + {C c AZ c )/r: 



(s k ,\V c (AZ c ,r)\s k ) dk' 



Sk> 



V c (AZ c ,r) 



C r AZ r 



s k ) dk' 



C c AZ C 



[fori 



[(k + k) ln(k max + k) 
k) ln(k max - k) - 2fcln(fc)] 



(18) 



AI(N GL , k is ,k max ) 
I GL (N GL ,h s ,k 

■max 

) - 1 (hi, , k max ) 



max \I(ki> , k max )\ 



)\0<i'<N GL 



(21) 



AI(N GL , 
and kn 



, ki s , k max ) is plotted for various values of Nql 
-"max on Fig. (HJ. One can see that it is very small, 
of the order of 1CT 6 to 10" 4 for N GL = 50, and de- 
creases by an order of magnitude when Nql is mul- 
tiplied by a factor of two. Interestingly, the shape of 
AI(N G L,h s , kmax) remains identical when N G l is mod- 
ified, with a and k max kept fixed. Precision is comparable 
to that of point-particle Coulomb potential for a = 0.45 
and 0.65 fin -1 , the values of interest for the lightest nu- 
clei, which are the most studied with the Berggren basis 
fl7H25j . As N GL rarely exceeds 100 in practical applica- 
tions and k max is typically of the order of 4 fm -1 [25|, the 
integral discretization inspired from Eq. (|16p is justified 
with Gauss-Legendre quadrature. 

We can now proceed to the discretization scheme we 
advocate for the general Berggren basis from a subtrac- 
tion method similar to that of Eq. (fT3"]) . The notations 
used will be identical to those of Eq. (fI3|) and Sec. Ill Bl 
We will also use the following definition: 

uf{r) 



N res < i < N, 



(22) 



where kf is defined in Eq. (|20|) . The subtraction scheme 
from which we start is the following: 

Cfe' (u k >\V c (AZ c ,r)\u k ) dk' 



L + 



(cfe/ Ufc'(r) u k (r) 



As we consider three finite a values for V c (AZ c ,r) [see 
Eq. 0], namely a = 0.25, 0.45 and 0.65 fm -1 , which 
correspond respectively to a number of nucleons equal 
to about 70, 10 and 5, calculating the radial integrals 
involving [V c (AZ c ,r) - (C c AZ c )/r] in Eq. JTSJ with 
300 Gauss-Legendre points defined on [0 : R], with R = 
30 fm, guarantees an almost exact reproduction of their 
numerical value. The Gauss-Legendre approximation of 
Eq. (|17l) then follows from the transformation mentioned 
above: 

Ngl-1 

lGL(N GL ,k lg ,k max ) = ^ (sk z ,jV c (AZ c ,r)\s kta ) Wi> s 

i 's=° 

+ (s k +\V c (AZ c ,r)\s k -) w ls , (19) 
where < i s < N G l and kf is defined as: 

kf = k is ± ^ , < i, < N GL . (20) 

47T 

In order to assess the precision of Eq. (|19[) when com- 
pared to Eq. (TT7)) . we define the following relative differ- 
ence: 



Cfc Sfc'(r) s fe (r)) V c (AZ Cl r) dr 

c fe / (s k ,\V c {AZ c ,r)\s k ) dk', 
Jl+ 



dk' 



(23) 



where Coulomb point potential is no longer utilized, con- 
trary to that of Eq. (fT5|i . Based on the latter theoretical 
and numerical arguments, the discretization scheme of 
Eq. (JTSJI becomes for Eq. (123)) : 



N-l 



(H>{ui'\V c {AZ c ,r)\ui) = E a , < i < N r . 



i'=0 



c, 



N-l 



+ / ( Ul (r) 2 - w ls s fela (r) 2 ) V C (AZ C , r) dr 



o 



s k+ {r)V c (AZ c ,r)s k - (r) dr 



Ci>(ui>\V c (AZ c ,r)\ui) =Ea, N res <i<N. (24) 



i'=0 



All off-diagonal matrix elements involving Sk is (r) 
and Sk'. (t) issued from the integrals of the 
right-hand side of Eq. ([25]) cancel out. Indeed, 



+00 



(cfe/ Ufe'(r) u fc (r) - Cfe s fc /(r) s k (r)) V c (AZ Cl r) dr 

can always be rewritten as c k > (u k >\V c (AZ c ,r)\u k ) — 
Cfe (s k i\V c (AZ c , r)\s k ), where the two latter matrix 
elements can be calculated with complex rotation of r 
as k 7^ k' (see Sec. Ill Ap . One can verify that replacing 

u % (r) 2 by uf(r)u~(r) and s kis (r) 2 by s k +(r)s k -(r) 

[see Eqs. (|2T)]) and (|22p] in the radial integral of the 
second line of Eq. (|2~4"|) results in an error of the order of 
u>f s /(47r) 2 , which can be neglected. Because kf ^ k^ , 
the cancellation we have noticed to take place for 
off-diagonal matrix elements now occurs in the modified 
diagonal radial integral. The off-diagonal discretization 
scheme then reads: 

AT-l 

c t e l + Ci> (ui>\V c (AZ c ,r)\ui) = E Ci , < i < N res 

i'=Q 

c t [ ei + (uf\V c (AZ c ,r)\Ui)] 

N-l 

+ a>(ui'\V c {AZ c ,r)\ui) =Ea, N res <i<N. (25) 

i'=0 

i' 

Note that the k max dependence of Eq. (fTSJ) has disap- 
peared in Eq. (f2"5j) . so that this scheme can be used in 
principle when k max — > +00. 
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FIG. 1: Relative difference AI(NaL,k,k ma x) [see Eq. (|21[) ] as a function of fc. A^gl is the number of Gauss-Legendre points 
in lGL(NcL,k,kmax) [see Eq. (fT§)) ] and fc varies from fc = to k = k max . We consider JVgl = 50, 100 and 200, k ma x = 1, 2 
and 4 fm -1 and a = 0.25, 0.45 and 0.65 fm -1 , where a defines the Coulomb potential V c (AZ c ,r) used [see Eq. Q]. Results 
are independent of AZ C as it simplifies in AI(NaL,k,k m ax) [see Eq. I|21|)]. Long-dashed, dashed, and dotted lines correspond 
respectively to a = 0.25, 0.45 and 0.65 fm -1 . Solid lines refer to the use of a point-particle Coulomb potential (C c AZ c )/r 
instead of V c (AZ c ,r) in Eqs. (|4]), (|17[) . and (|19[) . Data corresponding to kmax = 4 fm - and a — 0.25, 0.45 and 0.65 fm -1 have 
been scaled down by a factor of five, as indicated by "x 0.2 (a)". 



III. NUMERICAL EXAMPLES 

We consider now the expansion of the three single- 
particle states of the sd-shell with a Berggren basis, 
namely proton lsi/2, 0d$/2 and Qd 3 / 2 wave functions. 
The parameters defining the Woods-Saxon potential 
V\vs{ r ) an d Coulomb potential V c (Z c ,r) [see Eqs. §5§ 
and ([I])], common for basis and diagonalized Hamiltoni- 
ans [see Eq. ©], are d = 0.65 fm, R = 3 fm, V a = 52 



MeV, V so = 5 McV and a = 3^F/(4i? )- The latter 
value for a is chosen in order for V C (Z C , r = 0) to sustain 
the same value as that of the Coulomb potential defined 
from a uniformly charged-sphere of radius Rq (2f| . Basis 
and diagonalized Hamiltonians differ through the value 
of Z c in Eq. Q, in which z[ b) = 10 and Z ( c d) = 8 (see 
Sec. Ill Al for notations). These parameters are typical of 
Woods-Saxon potentials mimicking nuclei of A ~ 20 nu- 
cleons. The energies and widths of proton states of the 
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basis and diagonalized Hamiltonians, calculated with di- 
rect integration, are given in Table HI Note that the pro- 
ton Osi/2 state appears in the Si/ 2 Berggren basis, but, 
as it is well bound, the proton Osi/2 state just bears per- 
turbativc modification from basis to diagonalized Hamil- 
tonian, so that it is not interesting to study its case. It 
is important to state that it is not our aim to repro- 
duce experimental data, but only to study how precise 
Berggren basis expansion is if an infinite-range Coulomb 
potential is present in the diagonalized residual interac- 
tion. When wielding the cut method (see Sec. Ill C[) . the 
cut radius R of Eq. will be fixed at R = 75 fm for 
Sx/2 and d 5 / 2 partial waves, and at R — 35 fm for the 
(I3/2 partial wave, which are the values yielding the best 
precision for the cut method, while the Berggren basis 
contour [see Eq. (JTJ] consists of three segments of the 
complex /c-plane, delimited by the four points k m i n = 
fin. , k = 0.25-0. li fm -1 (siy 2 and d 5/ / 2 partial waves) 
or k = 0.4-0.39i fm -1 ((I3/2 partial wave), k = 1 fm -1 , 
and k ma x = 4 fm -1 . 



of their real and imaginary parts, defined by: 



The radial contour employed for integration with com- 
plex rotation of r in Eq. (JTJ is defined by fixing R — 
15 fm and 6 = —3tt/4, — 7r/4, tt/A or 37r/4. Improper 
integrals of Eq. (|7J| are indeed guaranteed to converge 
employing one of these angular values. The Berggren 
basis contour used in Eq. ([T]) is very similar to that 
of the cut method, except that k m i n is chosen so that 
\Ftn(kmi n R)\ + \k min F^(k min R)\ = 10~ 5 . Indeed, it is 
necessary for k m i n to be stricly positive, because, on the 
one hand, proton scattering states are very close to reg- 
ular Coulomb wave functions y/2/wFi v (kr) when k — > 
very small for moderate values of r, and on the other 
hand, H^(kr) functions enter integration effected with 
complex rotation of r in Eq. ([?])■ which bear very large 
modulus for moderate values of r. As a consequence, very 
important numerical cancellations occur between the dif- 
ferent improper integrals of Eq. 0. Nevertheless, proton 
scattering states with very small linear momentum play 
virtually no role in Berggren basis completeness, because 
of the extreme smallness of their amplitudes. Hence, pro- 
ton scattering states become important only when fc m j„ 
become sufficiently large, and the condition above has 
been shown to mitigate the numerical instability occur- 
ring for k ~ while yielding precise results. 



Comparison of energies and widths of considered pro- 
ton states provided by diagonalization to those deduced 
from direct integration are depicted in Tables UH IIII1 and 
IIV1 Precision of radial wave functions has also been 
taken into account by calculating the root mean squares 



»(«[«]) = \ 



N 



(»[«(r<)] - »[«e(ri)]) 5 



N 



,($*[«]) = ^ 



N 



(26) 



(27) 



where N is taken equal to 512, n = i ■ (R/N), for 
1 < i < N, is a set of uniformly distributed radii of 
[0 : R], u(r) is the diagonalized wave function, and u e (r) 
is the exact wave function, issued from direct integration. 
Root mean squares are illustrated in Fig. |2] for all studied 
cases, (i.e. the cut method, the subtraction method and 
the off-diagonal method; see Sees. IIIC| IIID| and III E|) . 
One can see that the cut method produces very poor re- 
sults, as it does not even reach the precision acquired 
with the off-diagonal method when the smallest value 
of the number of Gauss-Legendre scattering states Nql 
(see Sec. Ill Bp is employed. Moreover, for the 0g? 3 / 2 pro- 
ton state, although a rather good description of energy 
and width occurs, the reproduction of the wave func- 
tion is mediocre. Added to that, the choice of the cut 
radius could only be effected by comparison with exact 
results, whereas it has been checked that the two other 
methods are very robust when the L + contour param- 
eters are changed. The best reproduction of the con- 
sidered proton states clearly arises with the off-diagonal 
method, for energies, widths and wave functions. The 
subtraction method, while not being completely inaccu- 
rate, saturates very quickly to a wrong value for energies, 
widths and wave functions, when Nql increases. On the 
contrary, exponential convergence occurs with the off- 
diagonal method for both real and imaginary parts of the 
wave function when Nql augments. This is an intrigu- 
ing phenomenon, as the subtraction scheme of Eq. (|15p 
is based on an exact calculation of the integral exhibiting 
singularities, leaving a well-defined function to be inte- 
grated numerically, whereas that of Eq. (|25[) , which could 
be expected at best to be comparable to the subtraction 
method (see Sec. lIIDI and lIIE[) . surprisingly surpasses the 
latter by several orders of magnitudes (see Fig. [5])). It can 
be explained by noticing that the integrand of the first 
integral of the right-hand side of Eq. (H"3l) . constituted 
by radial integrals converging with complex rotation of 
r [see Eq. ([7)1] ■ while everywhere finite, is not analytic 
at kl ~ k. Indeed, it is equivalent, up to an unimpor- 
tant constant, to (k — k') ln(fc — k'), which does not even 
possess a finite derivative with respect to k' at k' = k. 
This implies that a Gauss-Legendre discretization of this 
integral will be far less precise than that of analytic func- 
tions, which can usually be well approximated by poly- 
nomials. On the contrary, the use of (uf\V c (AZ c , r)\u~) 
in Eq. ([23]) effectively replaces in Eq. © the function 
of k' (uk>\V c (AZ c ,r)\uk) , singular at kl = k, by an an- 
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TABLE I: Energies and widths of proton lsi/2, Ods/2 and CW3/2 states for both basis and diagonalized Hamiltonians, determined 
with direct integration. Energies are given in MeV and widths in keV. The basis states are denoted as "basis" and the 
diagonalized states as "diag." 





E basis (MeV) 


T basis (keV) 


E diag (MeV) 


T diag (keV) 


lSl/2 


1.09747 


134.623 


0.463324 


8.96828 


0d 5 /2 


1.48359 


11.9527 


0.666208 


0.525611 


0d 3 /2 


5.07435 


1353.51 


4.3003 


1091.3 



TABLE II: Energies and widths of proton lsi/2 state issued from the diagonalization of the Hamiltonian of Eq. ([2j with Z c 
= 8, expanded from a Berggren basis generated by a Hamiltonian of the same structure but bearing zi^ = 10 [see Eqs. (O 
and Q and Sec. Mil for the values of the other Hamiltonian parameters], as a function of the number of scattering states Ngl 
of the L + contour of Eq. {TJ, discretized with Gauss-Legendre quadrature. E designates the energy of the proton state, given 
in MeV, and F the width of the proton state, given in keV, while "exact" refers to results obtained with direct integration. 
The different discretization methods [i.e. the cut method, the subtraction method and the off-diagonal method (see Sees. Ill Cl 
III Dl and III E| )]. are denoted respectively as "cut," "sub," and "off-diag." 



Ngl 


E cut (MeV) 


F cut (keV) 


E sub (MeV) 


T sub (keV) 


E off-diag (MeV) 


T off-diag (keV) 


15 


0.461875 


-11.6596 


0.464574 


9.19011 


0.46396 


10.2211 


30 


0.465707 


13.4833 


0.463777 


8.26812 


0.463343 


8.97219 


45 


0.463476 


8.71097 


0.463709 


8.33267 


0.463334 


8.96171 


60 


0.463307 


8.68396 


0.463681 


8.36454 


0.463329 


8.96458 


75 


0.463227 


8.70558 


0.463667 


8.38006 


0.463328 


8.96595 


90 


0.46284 


8.88896 


0.463659 


8.3888 


0.463327 


8.96669 


105 


0.462952 


8.69106 


0.463654 


8.39421 


0.463326 


8.96712 


120 


0.462949 


8.62468 


0.46365 


8.3978 


0.463326 


8.9674 


exact 


0.463324 


8.96828 


0.463324 


8.96828 


0.463324 


8.96828 



alytic function. Gauss-Legendre quadrature then yields 
fast convergence to the numerical value of the integral 
defined by the off-diagonal method, which happens to be 
almost equal to the exact singular kernel (see Sec. Ill El) . 



IV. CONCLUSION 

Bases carrying a continuous part, such as Berggren 
bases, are very interesting as they allow the expansion 
of complex many-body wave functions of loosely bound 
and unbound states, as both proper asymptotic proper- 
ties and particle intercorrelations via configuration mix- 
ing are present therein. However, the use of Berggren 
bases is accompanied by mathematical difficulties gen- 
erated by the unbound character of the considered one- 
body states. Indeed, their completeness properties re- 
quest many more efforts to be proved than for discrete 
sets of states, to which the whole apparatus of the theory 
of compact operators can be applied. Furthermore, the 
very question of their numerical implementation is chal- 
lenging in the context of charged particles, as Coulomb 
wave functions have to be calculated for this type of 
problem. Discretization of Berggren bases is also central 
in numerical applications, which is conveniently effected 
with Gauss quadrature but demands discretization error 



to be assessed. Handling of infinite-range interactions is 
therefore problematic with unbound bases, because the 
representation of the former give rise to operators bearing 
singularities, which introduce infinities when Berggren 
bases are discretized. 

In the case of Coulomb potential, kernels are singular 
but integrable, so that frameworks built from subtraction 
techniques, based on analytical integration of singulari- 
ties, have been devised in the context of Fourier-Bessel 
transform. However, they could not be directly applied to 
Berggren bases, because they rely on the analytical char- 
acter of the Fourier-Bessel transform of Coulomb point 
potential. Hence, in this paper, three different discretiza- 
tion schemes have been studied when the Coulomb poten- 
tial is included in the Hamiltonian to diagonalize, namely 
the cut method, where the Coulomb potential is sup- 
pressed after a finite radius R, the subtraction method, 
similar to those utilized with Fourier-Bessel transform, 
and a new framework, the off-diagonal method, which 
amounts to substituting diagonal infinite matrix elements 
by close but finite off-diagonal matrix elements. 

Numerical applications have been considered for three 
different partial waves, with the examples of lsi/2, 0d 5 /2 
and 0d s / 2 resonant proton states, arising typically from 
studies of nuclei having A ~ 20 nucleons. Woods- 
Saxon potentials carrying different charges have been em- 
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TABLE III: Same as Table |TT1 but for the proton 0d 5 / 2 state. 



Ngl 


E cut (MeV) 


T cut (keV) 


E sub (MeV) 


T sub (keV) 


E off-diag (MeV) 


r off-diag (keV) 


15 


0.664431 


3.68108 


0.666482 


0.234635 


0.666428 


0.174831 


30 


0.665162 


-6.30699 


0.666251 


0.586056 


0.666216 


0.525826 


45 


0.66623 


0.502731 


0.66624 


0.578586 


0.666209 


0.527655 


60 


0.666207 


0.52976 


0.666237 


0.573284 


0.666209 


0.526746 


75 


0.666209 


0.536367 


0.666236 


0.5705 


0.666209 


0.526332 


90 


0.666206 


0.484164 


0.666235 


0.568849 


0.666209 


0.52611 


105 


0.666225 


0.502034 


0.666234 


0.567786 


0.666209 


0.525978 


120 


0.666226 


0.539324 


0.666234 


0.567058 


0.666209 


0.525892 


exact 


0.666208 


0.525611 


0.666208 


0.525611 


0.666208 


0.525611 



TABLE IV: Same as Table [TU but for the proton 0d 3/2 state. 



Ngl 


E cut (MeV) 


T cut (keV) 


E sub (MeV) 


r sub (keV) 


E off-diag (MeV) 


L off-diag (keV) 


15 


4.29692 


1091.1 


4.30016 


1091.48 


4.30017 


1091.49 


30 


4.30301 


1082.24 


4.3003 


1091.29 


4.30031 


1091.3 


45 


4.30016 


1091.42 


4.3003 


1091.29 


4.3003 


1091.3 


60 


4.30064 


1091.92 


4.3003 


1091.3 


4.3003 


1091.3 


75 


4.30072 


1092.15 


4.3003 


1091.3 


4.3003 


1091.3 


90 


4.30064 


1091.93 


4.3003 


1091.3 


4.3003 


1091.3 


105 


4.30069 


1091.41 


4.3003 


1091.3 


4.3003 


1091.3 


120 


4.30053 


1091.58 


4.3003 


1091.3 


4.3003 


1091.3 


exact 


4.3003 


1091.3 


4.3003 


1091.3 


4.3003 


1091.3 



ployed for basis-generating and diagonalized Hamiltoni- 
ans. While the cut method could be expected to convey 
poor precision, the off-diagonal method has been demon- 
strated to outperform the standard subtraction method. 
This has been explained by the fact that Gauss-Legendre 
quadrature is applied within the subtraction method 
to integrate finite but nonanalytic functions, whereas 
only smooth integrands arc treated with Gauss-Legendre 
quadrature within the off-diagonal method. According to 



the present study, the latter technique should be taken 
into account seriously when diagonalizing Hamiltonians 
possessing an infinite-range Coulomb part. 
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